Análisis de clases latentes exploratoria y comparativa sin y con predictores, sin pueblo originario, sin año y recategorización de edad mujer y semana gestacional hito 1, caracterización de clases y medidas de ajuste
script src = "https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js" $(document).ready(function() {
$('body').prepend('<div class=\"zoomDiv\"><img src=\"\" class=\"zoomImg\"></div>');
// onClick function for all plots (img's)
$('img:not(.zoomImg)').click(function() {
$('.zoomImg').attr('src', $(this).attr('src')).css({width: '100%'});
$('.zoomDiv').css({opacity: '1', width: 'auto', border: '1px solid white', borderRadius: '5px', position: 'fixed', top: '50%', left: '50%', marginRight: '-50%', transform: 'translate(-50%, -50%)', boxShadow: '0px 0px 50px #888888', zIndex: '50', overflow: 'auto', maxHeight: '100%'});
});
// onClick function for zoomImg
$('img.zoomImg').click(function() {
$('.zoomDiv').css({opacity: '0', width: '0%'});
});
});
<script src="hideOutput.js"></script> $(document).ready(function() {
$chunks = $('.fold');
$chunks.each(function () { // add button to source code chunks
if ( $(this).hasClass('s') ) {
$('pre.r', this).prepend("<div class=\"showopt\">Show Source</div><br style=\"line-height:22px;\"/>");
$('pre.r', this).children('code').attr('class', 'folded');
} // add button to output chunks
if ( $(this).hasClass('o') ) {
$('pre:not(.r)', this).has('code').prepend("<div class=\"showopt\">Show Output</div><br style=\"line-height:22px;\"/>");
$('pre:not(.r)', this).children('code:not(r)').addClass('folded'); // add button to plots
$(this).find('img').wrap('<pre class=\"plot\"></pre>');
$('pre.plot', this).prepend("<div class=\"showopt\">Show Plot</div><br style=\"line-height:22px;\"/>");
$('pre.plot', this).children('img').addClass('folded');
}
}); // hide all chunks when document is loaded
$('.folded').css('display', 'none') // function to toggle the visibility
$('.showopt').click(function() {
var label = $(this).html();
if (label.indexOf("Show") >= 0) {
$(this).html(label.replace("Show", "Hide"));
} else {
$(this).html(label.replace("Hide", "Show"));
}
$(this).siblings('code, img').slideToggle('fast', 'swing');
});
}); Cargamos los datos
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 653061 34.9 1307962 69.9 929440 49.7
Vcells 1172897 9.0 8388608 64.0 2056718 15.7
load("data2_lca2_adj_sin_po_ano.RData")
Cargamos los paquetes
knitr::opts_chunk$set(echo = TRUE)
if(!require(poLCA)){install.packages("poLCA")}
if(!require(poLCAParallel)){devtools::install_github("QMUL/poLCAParallel@package")}
if(!require(compareGroups)){install.packages("compareGroups")}
if(!require(parallel)){install.packages("parallel")}
if(!require(Hmisc)){install.packages("Hmisc")}
if(!require(tidyverse)){install.packages("tidyverse")}
try(if(!require(sjPlot)){install.packages("sjPlot")})
if(!require(emmeans)){install.packages("emmeans")}
if(!require(nnet)){install.packages("nnet")}
if(!require(here)){install.packages("here")}
if(!require(doParallel)){install.packages("doParallel")}
if(!require(progress)){install.packages("progress")}
if(!require(caret)){install.packages("caret")}
if(!require(rpart)){install.packages("rpart")}
if(!require(rpart.plot)){install.packages("rpart.plot")}
if(!require(partykit)){install.packages("partykit")}
if(!require(randomForest)){install.packages("randomForest")}
if(!require(ggcorrplot)){install.packages("ggcorrplot")}
if(!require(polycor)){install.packages("polycor")}
if(!require(tableone)){install.packages("tableone")}
if(!require(broom)){install.packages("broom")}
if(!require(plotly)){install.packages("plotly")}
if(!require(rsvg)){install.packages("rsvg")}
if(!require(DiagrammeRsvg)){install.packages("DiagrammeRsvg")}
if(!require(effects)){install.packages("effects")}
#Si probs.start se establece en NULL (predeterminado) al llamar Polca, a continuación, la función genera los valores de partida al azar en cada ejecución. Esto significa que repite carreras de Polca normalmente producirán resultados con los mismos parámetros estimados (correspondiente a la misma el máximo diario de probabilidad), pero con etiquetas de clase latentes reordenados
#https://drive.google.com/file/d/10njMh5JEcqaBgCnoZdJ1uk3uCEkocDez/view?usp=share_link
#http://daob.nl/wp-content/uploads/2015/07/ESRA-course-slides.pdf
#https://docs.google.com/document/d/1LVeDpAP6CfT3n8B6HhHcc_SRnzO0JBoT/edit
#bvr(ppio_newList$lc_entropy_best_model)
#A list of matrices of class-conditional response probabilities to be used as the starting values for the estimation algorithm. Each matrix in the list corresponds to one manifest variable, with one row for each latent class, and one column for each outcome. The default is NULL, producing random starting values. Note that if nrep>1, then any user-specified probs.start values are only used in the first of the nrep attempts.
#The poLCA.reorder function takes as its first argument the list of starting values probs.start, and as its second argument a vector describing the desired reordering of the latent classes.
new.probs.start <- poLCA.reorder(LCA_best_model_ppio$probs.start, order(LCA_best_model_ppio$P, decreasing = TRUE))
#new.probs.start <-poLCA.reorder(probs.start,c(4,1,3,2))
#A continuación, ejecute PoLCA, una vez más, esta vez utilizando los valores iniciales reordenados en la llamada de función.
#The argument nrep=5 tells the program to repeat nrep times, and keep the fit with the highest likelihood to try to avoid local maxima.
#.maxiter – The maximum number of iterations through which the estimation algorithm will cycle.
#.nrep - Number of times to estimate the model, using different values of probs.start. (default is one)
set.seed(2125)
LCA_best_model_mod<-
poLCA(f_preds, mydata_preds3, nclass=length(LCA_best_model_ppio$P),
maxiter=10000,tol=1e-5, na.rm=FALSE,nrep=1e3, verbose=TRUE, calc.se=TRUE,probs.start=new.probs.start)
Conditional item response (column) probabilities,
by outcome variable, for each class (row)
$CAUSAL
Pr(1) Pr(2) Pr(3) Pr(4)
class 1: 0 0.0233 0.0000 0.9767
class 2: 0 0.1928 0.8072 0.0000
class 3: 0 0.4647 0.5353 0.0000
class 4: 0 0.4038 0.5962 0.0000
class 5: 0 0.0098 0.0000 0.9902
$EDAD_MUJER_REC
Pr(1) Pr(2) Pr(3) Pr(4) Pr(5)
class 1: 0.0000 0.1649 0.2335 0.4515 0.1501
class 2: 0.0000 0.0000 0.0405 0.4255 0.5340
class 3: 0.0043 0.0085 0.1703 0.5025 0.3144
class 4: 0.0072 0.0180 0.2106 0.4571 0.3071
class 5: 0.0017 0.3472 0.2463 0.2958 0.1089
$PAIS_ORIGEN_REC
Pr(1) Pr(2) Pr(3)
class 1: 0.0000 0.0001 0.9999
class 2: 0.0000 0.9252 0.0748
class 3: 0.0339 0.0431 0.9230
class 4: 0.0002 0.9998 0.0000
class 5: 0.0000 0.9935 0.0065
$HITO1_EDAD_GEST_SEM_REC
Pr(1) Pr(2) Pr(3) Pr(4)
class 1: 0.0191 0.9809 0.0000 0.0000
class 2: 0.0150 0.3653 0.5698 0.0499
class 3: 0.0230 0.1163 0.7339 0.1268
class 4: 0.0291 0.1767 0.6594 0.1348
class 5: 0.0089 0.9877 0.0035 0.0000
$MACROZONA
Pr(1) Pr(2) Pr(3) Pr(4) Pr(5) Pr(6)
class 1: 0.0062 0.5562 0.0870 0.0166 0.2960 0.0379
class 2: 0.0000 0.7147 0.0922 0.0556 0.0698 0.0677
class 3: 0.0059 0.5683 0.1051 0.0594 0.2083 0.0530
class 4: 0.0022 0.3208 0.1932 0.2004 0.0935 0.1899
class 5: 0.0042 0.3849 0.1747 0.1452 0.0862 0.2048
$PREV_TRAMO_REC
Pr(1) Pr(2) Pr(3) Pr(4) Pr(5)
class 1: 0.0253 0.0062 0.5225 0.1669 0.2791
class 2: 0.0000 0.8096 0.0000 0.1839 0.0065
class 3: 0.0123 0.0139 0.6086 0.2871 0.0780
class 4: 0.0012 0.0251 0.6227 0.3474 0.0036
class 5: 0.0017 0.0695 0.7236 0.1998 0.0053
Estimated class population shares
0.0425 0.1265 0.1371 0.541 0.1529
Predicted class memberships (by modal posterior prob.)
0.0425 0.118 0.1325 0.5563 0.1507
=========================================================
Fit for 5 latent classes:
=========================================================
number of observations: 3789
number of estimated parameters: 109
residual degrees of freedom: 3680
maximum log-likelihood: -22209.52
AIC(5): 44637.03
BIC(5): 45317.18
G^2(5): 1632.458 (Likelihood ratio/deviance statistic)
X^2(5): 695543.1 (Chi-square goodness of fit)
output_LCA_best_model_mod<-capture.output(LCA_best_model_mod)
glance_LCA_best_model_mod<-glance(LCA_best_model_mod)
mydata_preds_LCA1 <- augment(LCA_best_model_mod, data = mydata_preds3)
# fig.height=15,
## If you are interested in the population-shares of the classes, you can get them like this:
warning(paste("Probabilidades posteriores: ",
paste(round(colMeans(LCA_best_model_mod$posterior)*100,2), collapse=", ")
))
## or you inspect the estimated class memberships:
warning(paste("Probabildiades predichas: ",
paste(round(prop.table(table(LCA_best_model_mod$predclass)),4)*100, collapse=", ")
))
traductor_cats <- readxl::read_excel("tabla12_corr.xlsx") %>%
dplyr::mutate(level=readr::parse_double(level)) %>%
#dplyr::mutate(level= dplyr::case_when(grepl("CAUSAL",Name)~ level-1,T~level)) %>%
#dplyr::mutate(level= dplyr::case_when(grepl("AÑO",Name)~ level-1,T~level)) %>%
dplyr::mutate(charactersitic=gsub(" \\(%\\)", "", Name))
lcmodel <- reshape2::melt(LCA_best_model_mod$probs, level=2)
lcmodel<- lcmodel %>%
dplyr::mutate(pr=as.numeric(gsub("[^0-9.]+", "", Var2))) %>%
dplyr::left_join(traductor_cats[,c("charactersitic", "level", "CATEGORIA")], by= c("L2"="charactersitic", "pr"="level"))
lcmodel$text_label<-paste0("Categoria:",lcmodel$CATEGORIA,"<br>%: ",scales::percent(lcmodel$value))
zp1 <- ggplot(lcmodel,aes(x = L2, y = value, fill = Var2, label=text_label))
zp1 <- zp1 + geom_bar(stat = "identity", position = "stack")
zp1 <- zp1 + facet_grid(Var1 ~ .)
zp1 <- zp1 + scale_fill_brewer(type="seq", palette="Greys", na.value = "white") +theme_bw()
zp1 <- zp1 + labs(y = "Porcentaje de probabilidad de respuesta",
x = "",
fill ="Cateorías de\nRespuesta")
zp1 <- zp1 + theme( axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
panel.grid.major.y=element_blank())
zp1 <- zp1 + guides(fill = guide_legend(reverse=TRUE))
zp1 <- zp1 + theme(axis.text.x = element_text(angle = 30, hjust = 1))
#print(zp1)
ggplotly(zp1, tooltip = c("text_label"))%>% layout(xaxis= list(showticklabels = T),height=600, width=800)
Figure 1: Selected Model
#In this case, residuals are actual cell counts vs. expected cell counts.
bvr_LCA_best_model_mod<-bvr(LCA_best_model_mod)
#conditional probabilities
#Pr(B1=1|Class 3)
posteriors <- data.frame(LCA_best_model_mod$posterior, predclass=LCA_best_model_mod$predclass)
classification_table <- plyr::ddply(posteriors, "predclass", function(x) colSums(x[,1:length(LCA_best_model_ppio$P)]))
clasification_errors<- 1 - sum(diag(as.matrix(classification_table[,2:(length(LCA_best_model_ppio$P)+1)]))) / sum(classification_table[,2:(length(LCA_best_model_ppio$P)+1)])
warning(paste("Error de clasificación: ", round(clasification_errors,2)))
Warning: Error de clasificación: 0.05
entropy_alt <- function(p) sum(-p * log(p))
error_prior <- entropy_alt(LCA_best_model_mod$P) # Class proportions
error_post <- mean(apply(LCA_best_model_mod$posterior, 1, entropy_alt),na.rm=T)
R2_entropy_alt <- (error_prior - error_post) / error_prior
warning(paste("Entropía: ", round(R2_entropy_alt,2)))
Warning: Entropía: 0.81
#https://stackoverflow.com/questions/72783185/entropy-calculation-gives-nan-is-applying-na-omit-a-valid-tweak
entropy.safe <- function (p) {
if (any(p > 1 | p < 0)) stop("probability must be between 0 and 1")
log.p <- numeric(length(p))
safe <- p != 0
log.p[safe] <- log(p[safe])
sum(-p * log.p)
}
error_prior2 <- entropy.safe(LCA_best_model_mod$P) # Class proportions
error_post2 <- mean(apply(LCA_best_model_mod$posterior, 1, entropy.safe),na.rm=T)
R2_entropy_safe <- (error_prior2 - error_post2) / error_prior2
warning(paste("Entropía segura: ", round(R2_entropy_safe,2)))
Warning: Entropía segura: 0.89
entropy.brutal <- function (p) {
if (any(p > 1 | p < 0)) stop("probability must be between 0 and 1")
log.p <- log(p)
## as same as sum(na.omit(-p * log.p))
sum(-p * log.p, na.rm = TRUE)
}
error_prior3 <- entropy.brutal(LCA_best_model_mod$P) # Class proportions
error_post3 <- mean(apply(LCA_best_model_mod$posterior, 1, entropy.brutal),na.rm=T)
R2_entropy_brutal <- (error_prior3 - error_post3) / error_prior3
warning(paste("Entropía brutal: ", round(R2_entropy_brutal,2)))
Warning: Entropía brutal: 0.89
Warning: Entropía (solución Oberski): 0.89
#\#minimum average posterior robability of cluster membership (\>0.7), interpretability (classes are clearly distinguishable), and parsimony (each class has a sufficient sample size for further analysis; n≥50).
Ver si la exclusión de casos que no calzan en alguna clase tiene consecuencias.
#To evaluate whether the exclusion of cases would bias the LCA results, a sensitivity analysis was carried out. We conducted T-Tests and Wilcoxon–Mann–Whitney tests (for non-parametric data) to compare included and excluded records in terms of demographic and clinical background characteristics and baseline pain scores (all 638 patients completed pain intensity, frequency and impairment scales).
tidy(LCA_best_model_mod) %>%
# dplyr::mutate(variable= dplyr::case_when(variable=="naq1"~"naq01",
# variable=="naq2"~"naq02",
# variable=="naq3"~"naq03",
# variable=="naq4"~"naq04",
# variable=="naq5"~"naq05",
# variable=="naq6"~"naq06",
# variable=="naq7"~"naq07",
# variable=="naq8"~"naq08",
# variable=="naq9"~"naq09",
# TRUE~variable)) %>%
ggplot(aes(outcome, estimate, color = factor(class), group = class)) +
geom_line() +
facet_wrap(~variable, nrow = 4) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
theme_bw()+
theme(legend.position = "bottom")
Figure 2: Modelo seleccionado
.
new.probs.start_adj <- poLCA.reorder(LCA_best_model_adj_ppio$probs.start,
order(LCA_best_model_adj_ppio$P, decreasing = TRUE))
set.seed(2125)
LCA_best_model_adj_mod<-
poLCA(f_adj, mydata_preds3, nclass=length(LCA_best_model_adj_ppio$P),
maxiter=10000,tol=1e-5, na.rm=FALSE,nrep=1e3, verbose=TRUE, calc.se=TRUE,probs.start=new.probs.start_adj)
Conditional item response (column) probabilities,
by outcome variable, for each class (row)
$CAUSAL
Pr(1) Pr(2) Pr(3) Pr(4)
class 1: 0 0.0049 0.0000 0.9951
class 2: 0 0.4626 0.5374 0.0000
class 3: 0 0.4068 0.5932 0.0000
class 4: 0 0.0247 0.0000 0.9753
class 5: 0 0.2021 0.7979 0.0000
$EDAD_MUJER_REC
Pr(1) Pr(2) Pr(3) Pr(4) Pr(5)
class 1: 0.0017 0.3483 0.2460 0.2951 0.1088
class 2: 0.0000 0.0089 0.1747 0.5022 0.3142
class 3: 0.0082 0.0181 0.2114 0.4563 0.3060
class 4: 0.0000 0.1648 0.2331 0.4511 0.1509
class 5: 0.0000 0.0000 0.0393 0.4346 0.5262
$PAIS_ORIGEN_REC
Pr(1) Pr(2) Pr(3)
class 1: 0.0000 0.9933 0.0067
class 2: 0.0121 0.0011 0.9868
class 3: 0.0059 0.9941 0.0000
class 4: 0.0000 0.0028 0.9972
class 5: 0.0000 0.9069 0.0931
$HITO1_EDAD_GEST_SEM_REC
Pr(1) Pr(2) Pr(3) Pr(4)
class 1: 0.0088 0.9877 0.0035 0.0000
class 2: 0.0204 0.1141 0.7361 0.1294
class 3: 0.0290 0.1733 0.6603 0.1374
class 4: 0.0192 0.9808 0.0000 0.0000
class 5: 0.0181 0.3703 0.5721 0.0396
$MACROZONA
Pr(1) Pr(2) Pr(3) Pr(4) Pr(5) Pr(6)
class 1: 0.0038 0.3850 0.1751 0.1453 0.0860 0.2048
class 2: 0.0000 0.5660 0.1060 0.0623 0.2147 0.0510
class 3: 0.0038 0.3203 0.1918 0.1999 0.0948 0.1894
class 4: 0.0062 0.5560 0.0871 0.0165 0.2957 0.0385
class 5: 0.0000 0.7205 0.0937 0.0497 0.0700 0.0660
$PREV_TRAMO_REC
Pr(1) Pr(2) Pr(3) Pr(4) Pr(5)
class 1: 0.0018 0.0697 0.7238 0.1996 0.0051
class 2: 0.0098 0.0007 0.6260 0.2819 0.0817
class 3: 0.0020 0.0305 0.6254 0.3384 0.0037
class 4: 0.0256 0.0061 0.5214 0.1674 0.2795
class 5: 0.0000 0.7578 0.0050 0.2282 0.0091
Estimated class population shares
0.1521 0.1251 0.5466 0.0426 0.1335
Predicted class memberships (by modal posterior prob.)
0.1507 0.1288 0.5585 0.0425 0.1196
=========================================================
Fit for 5 latent classes:
=========================================================
2 / 1
Coefficient Std. error t value Pr(>|t|)
(Intercept) 0.41261 0.23412 1.762 0.078
outcome1 -0.68583 0.21270 -3.224 0.001
=========================================================
3 / 1
Coefficient Std. error t value Pr(>|t|)
(Intercept) 2.25089 0.16625 13.540 0
outcome1 -1.13159 0.17346 -6.524 0
=========================================================
4 / 1
Coefficient Std. error t value Pr(>|t|)
(Intercept) -1.57566 0.40369 -3.903 0.000
outcome1 0.32738 0.40273 0.813 0.416
=========================================================
5 / 1
Coefficient Std. error t value Pr(>|t|)
(Intercept) -0.66059 0.37749 -1.750 0.080
outcome1 0.56682 0.34648 1.636 0.102
=========================================================
number of observations: 3789
number of estimated parameters: 113
residual degrees of freedom: 3676
maximum log-likelihood: -22149.64
AIC(5): 44525.27
BIC(5): 45230.38
X^2(5): 4907035 (Chi-square goodness of fit)
ALERT: estimation algorithm automatically restarted with new initial values
output_LCA_best_model_adj_mod<-capture.output(LCA_best_model_adj_mod)
glance_LCA_best_model_adj_mod<-glance(LCA_best_model_adj_mod)
mydata_preds_LCA2 <- augment(LCA_best_model_adj_mod, data = mydata_preds3)
# fig.height=15,
## If you are interested in the population-shares of the classes, you can get them like this:
warning(paste("Probabilidades posteriores: ",
paste(round(colMeans(LCA_best_model_adj_mod$posterior)*100,2), collapse=", ")
))
## or you inspect the estimated class memberships:
warning(paste("Probabildiades predichas: ",
paste(round(prop.table(table(LCA_best_model_adj_mod$predclass)),4)*100, collapse=", ")
))
lcmodel_adj <- reshape2::melt(LCA_best_model_adj_mod$probs, level=2)
lcmodel_adj<- lcmodel_adj %>%
dplyr::mutate(pr=as.numeric(gsub("[^0-9.]+", "", Var2))) %>%
dplyr::left_join(traductor_cats[,c("charactersitic", "level", "CATEGORIA")], by= c("L2"="charactersitic", "pr"="level"))
lcmodel_adj$text_label<-paste0("Categoria:",lcmodel_adj$CATEGORIA,"<br>%: ",scales::percent(lcmodel_adj$value))
zp2 <- ggplot(lcmodel_adj,aes(x = L2, y = value, fill = Var2, label=text_label))
zp2 <- zp2 + geom_bar(stat = "identity", position = "stack")
zp2 <- zp2 + facet_grid(Var1 ~ .)
zp2 <- zp2 + scale_fill_brewer(type="seq", palette="Greys", na.value = "white") +theme_bw()
zp2 <- zp2 + labs(y = "Porcentaje de probabilidad de respuesta",
x = "",
fill ="Cateorías de\nRespuesta")
zp2 <- zp2 + theme( axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
panel.grid.major.y=element_blank())
zp2 <- zp2 + guides(fill = guide_legend(reverse=TRUE))
zp2 <- zp2 + theme(axis.text.x = element_text(angle = 30, hjust = 1))
#print(zp1)
ggplotly(zp2, tooltip = c("text_label"))%>% layout(xaxis= list(showticklabels = T),height=600, width=800)
Figure 3: Selected Model
#In this case, residuals are actual cell counts vs. expected cell counts.
bvr_LCA_best_model_adj_mod<-bvr(LCA_best_model_adj_mod)
#conditional probabilities
#Pr(B1=1|Class 3)
posteriors_adj <- data.frame(LCA_best_model_adj_mod$posterior,
predclass=LCA_best_model_adj_mod$predclass)
classification_table_adj <- plyr::ddply(posteriors, "predclass", function(x) colSums(x[,1:length(LCA_best_model_adj_mod$P)]))
clasification_errors_adj<- 1 - sum(diag(as.matrix(classification_table_adj[,2:(length(LCA_best_model_adj_mod$P)+1)]))) / sum(classification_table_adj[,2:(length(LCA_best_model_adj_mod$P)+1)])
warning(paste("Error de clasificación: ", round(clasification_errors_adj,2)))
Warning: Error de clasificación: 0.05
entropy_alt <- function(p) sum(-p * log(p))
error_prior_adj <- entropy_alt(LCA_best_model_adj_mod$P) # Class proportions
error_post_adj <- mean(apply(LCA_best_model_adj_mod$posterior, 1, entropy_alt),na.rm=T)
R2_entropy_alt_adj <- (error_prior_adj - error_post_adj) / error_prior_adj
warning(paste("Entropía: ", round(R2_entropy_alt_adj,2)))
Warning: Entropía: 0.82
#https://stackoverflow.com/questions/72783185/entropy-calculation-gives-nan-is-applying-na-omit-a-valid-tweak
entropy.safe <- function (p) {
if (any(p > 1 | p < 0)) stop("probability must be between 0 and 1")
log.p <- numeric(length(p))
safe <- p != 0
log.p[safe] <- log(p[safe])
sum(-p * log.p)
}
error_prior2_adj <- entropy.safe(LCA_best_model_adj_mod$P) # Class proportions
error_post2_adj <- mean(apply(LCA_best_model_adj_mod$posterior, 1, entropy.safe),na.rm=T)
R2_entropy_safe_adj <- (error_prior2_adj - error_post2_adj) / error_prior2_adj
warning(paste("Entropía segura: ", round(R2_entropy_safe,2)))
Warning: Entropía segura: 0.89
Warning: Entropía (solución Oberski): 0.9
#\#minimum average posterior robability of cluster membership (\>0.7), interpretability (classes are clearly distinguishable), and parsimony (each class has a sufficient sample size for further analysis; n≥50).
Ver si la exclusión de casos que no calzan en alguna clase tiene consecuencias.
#To evaluate whether the exclusion of cases would bias the LCA results, a sensitivity analysis was carried out. We conducted T-Tests and Wilcoxon–Mann–Whitney tests (for non-parametric data) to compare included and excluded records in terms of demographic and clinical background characteristics and baseline pain scores (all 638 patients completed pain intensity, frequency and impairment scales).
tidy(LCA_best_model_adj_mod) %>%
# dplyr::mutate(variable= dplyr::case_when(variable=="naq1"~"naq01",
# variable=="naq2"~"naq02",
# variable=="naq3"~"naq03",
# variable=="naq4"~"naq04",
# variable=="naq5"~"naq05",
# variable=="naq6"~"naq06",
# variable=="naq7"~"naq07",
# variable=="naq8"~"naq08",
# variable=="naq9"~"naq09",
# TRUE~variable)) %>%
ggplot(aes(outcome, estimate, color = factor(class), group = class)) +
geom_line() +
facet_wrap(~variable, nrow = 4) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
theme_bw()+
theme(legend.position = "bottom")
Figure 4: Modelo seleccionado
# fig.cap= "Probabilidad predicha de presentar interrupción del embarazo por clase",
cbind(
effects::allEffects(LCA_best_model_adj_mod, se=T, confint=T)$outcome$prob, effects::allEffects(LCA_best_model_adj_mod, se=T, confint=T)$outcome$lower.prob,
effects::allEffects(LCA_best_model_adj_mod, se=T, confint=T)$outcome$upper.prob) %>%
data.frame() %>%
dplyr::mutate_all(~round(.,2))->df_model_LCA
# Define the number of probabilities
num_cols <- length(LCA_best_model_adj_mod$P.se)
df_model_LCA2<- data.frame(0)
for(i in 1:num_cols){
df_model_LCA2[,i] <- i
}
for (i in 1:num_cols) {
col_name <- paste0("prob.X", i)
L_name <- paste0("L.prob.X", i)
U_name <- paste0("U.prob.X", i)
new_col_name <- paste0("prob_c", i)
df_model_LCA2 <- cbind.data.frame(df_model_LCA2,df_model_LCA %>%
mutate(!!new_col_name := paste0(!!sym(col_name), "(95%CI=", !!sym(L_name), ",", !!sym(U_name),")")))
}
df_model_LCA2 %>% select(starts_with("prob_c")) %>%
t() %>% data.table::data.table(keep.rownames=T) %>%
dplyr::mutate(rn=gsub("prob_c","",rn)) %>%
knitr::kable("markdown", caption="Probabilidad de pertenecer a una clase, según interrupción del embarazo", col.names = c("Clase","No interrumpe", "Interrumpe"))
| Clase | No interrumpe | Interrumpe |
|---|---|---|
| 1 | 0.19(95%CI=0.16,0.23) | 0.18(95%CI=0.17,0.19) |
| 2 | 0.18(95%CI=0.15,0.21) | 0.15(95%CI=0.14,0.17) |
| 3 | 0.3(95%CI=0.26,0.33) | 0.39(95%CI=0.37,0.4) |
| 4 | 0.16(95%CI=0.13,0.19) | 0.13(95%CI=0.12,0.14) |
| 5 | 0.17(95%CI=0.15,0.21) | 0.15(95%CI=0.14,0.17) |
#########################################################
#########################################################
### Posterior probability calculation ###
### Assign class based on maximum probability ###
### Note: additional prep for Table1 package ###
### 1) Convert all categorical variables to ###
### factors ###
### 2) Continuous variables as numeric ###
### 3) Pull out number from census strings ###
#########################################################
#########################################################
save.image("data2_lca3_sin_po_ano.RData")
require(tidyverse)
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
) %>%
DT::datatable(filter = 'top', colnames = c('Row number' =1,'Variable' = 2, 'Percentage'= 3),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left;',
'', htmltools::em('Packages')),
options=list(
initComplete = htmlwidgets::JS(
"function(settings, json) {",
"$(this.api().tables().body()).css({
'font-family': 'Helvetica Neue',
'font-size': '50%',
'code-inline-font-size': '15%',
'white-space': 'nowrap',
'line-height': '0.75em',
'min-height': '0.5em'
});",#;
"}")))